#ifndef _BLIP2_H_
#define _BLIP2_H_

#include "Equations/TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class BLIP2;

  template <typename T>
  struct Typeinfo <BLIP2<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <BLIP2<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <BLIP2<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };

  
  template<typename Tprec, int Dim>
  class BLIP2 : public TwoPhaseEquation<BLIP2<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<BLIP2<Tprec, Dim> > TP_BLIP2;
    
    using GeneralEquation< TP_BLIP2 >::aE;
    using GeneralEquation< TP_BLIP2 >::aW;
    using GeneralEquation< TP_BLIP2 >::aN;
    using GeneralEquation< TP_BLIP2 >::aS;
    using GeneralEquation< TP_BLIP2 >::aB;
    using GeneralEquation< TP_BLIP2 >::aP;
    using GeneralEquation< TP_BLIP2 >::sp;
    using GeneralEquation< TP_BLIP2 >::dx;
    using GeneralEquation< TP_BLIP2 >::dy;
    using GeneralEquation< TP_BLIP2 >::dz;
    using GeneralEquation< TP_BLIP2 >::bi;
    using GeneralEquation< TP_BLIP2 >::ei;
    using GeneralEquation< TP_BLIP2 >::bj;
    using GeneralEquation< TP_BLIP2 >::ej;
    using GeneralEquation< TP_BLIP2 >::bk;
    using GeneralEquation< TP_BLIP2 >::ek;
    using GeneralEquation< TP_BLIP2 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_BLIP2 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_BLIP2 >::applyBoundaryConditions3D;

    using TP_BLIP2::S;
    using TP_BLIP2::phi_0;
    using TP_BLIP2::Srw;
    using TP_BLIP2::Sro;
    using TP_BLIP2::mu_w;
    using TP_BLIP2::mu_o;
    using TP_BLIP2::k;

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    BLIP2() : TwoPhaseEquation<BLIP2<prec_t, Dim > >() { }    
    ~BLIP2() { };
    
    inline bool calcCoefficients1D(); 
    inline bool calcCoefficients2D();
    inline bool calcCoefficients3D() { };
    inline void printInfo() { std::cout << " BLIP2 "; }
  };

/*
 *  Quadratic for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool BLIP2<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;

    // Quadratic
    static prec_t mult_o = k / 
      ( (1 - Srw - Sro) * (1 - Srw - Sro) * mu_o * dx );
    static prec_t mult_w = k / 
      ( (1 - Srw - Sro) * (1 - Srw - Sro) * mu_w * dx );
    
    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {      
      // Upwind
      if ( phi_0(i+1) >= phi_0(i) ) Sw_e = S(i+1);
      else                  Sw_e = S(i);
      if ( phi_0(i-1) >= phi_0(i) ) Sw_w = S(i-1);
      else                  Sw_w = S(i);
      
      // Quadratic
      aE (i) = (1 - Sro - Sw_e) * (1 - Sro - Sw_e) * mult_o 
	+ (Sw_e - Srw) * (Sw_e - Srw) * mult_w ;
      aW (i) = (1 - Sro - Sw_w) * (1 - Sro - Sw_w) * mult_o 
	+ (Sw_w - Srw) * (Sw_w - Srw) * mult_w ;

      aP (i) = aE (i) + aW (i);      
    }
    applyBoundaryConditions1D();
    return 0;
}

/*
 *  Quadratic for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool BLIP2<Tprec, Dim>::calcCoefficients2D () 
{
  static prec_t Sw_e, Sw_w,  Sw_n, Sw_s;

    // Quadratic
    static prec_t mult_o = k / 
      ( (1 - Srw - Sro) * (1 - Srw - Sro) * mu_o );
    static prec_t mult_w = k / 
      ( (1 - Srw - Sro) * (1 - Srw - Sro) * mu_w );

    static prec_t dx_dy = dx / dy;
    static prec_t dy_dx = dy / dx;   
    
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i)
      for (int j = bj; j <= ej; ++j) {  
	
	// Upwind
      if ( phi_0(i+1, j) >= phi_0(i, j) ) Sw_e = S(i+1, j);
      else                                Sw_e = S(i,   j);
      if ( phi_0(i-1, j) >= phi_0(i, j) ) Sw_w = S(i-1, j);
      else                                Sw_w = S(i,   j);

      if ( phi_0(i, j+1) >= phi_0(i, j) ) Sw_n = S(i, j+1);
      else                                Sw_n = S(i,   j);
      if ( phi_0(i, j-1) >= phi_0(i, j) ) Sw_s = S(i, j-1);
      else                                Sw_s = S(i,   j);

      
      // Quadratic
      aE (i,j) = ( (1 - Sro - Sw_e) * (1 - Sro - Sw_e) * mult_o 
		   + (Sw_e - Srw) * (Sw_e - Srw) * mult_w ) * dy_dx;
      aW (i,j) = ( (1 - Sro - Sw_w) * (1 - Sro - Sw_w) * mult_o 
		   + (Sw_w - Srw) * (Sw_w - Srw) * mult_w ) * dy_dx;
      aN (i,j) = ( (1 - Sro - Sw_n) * (1 - Sro - Sw_n) * mult_o 
		   + (Sw_n - Srw) * (Sw_n - Srw) * mult_w ) * dx_dy;
      aS (i,j) = ( (1 - Sro - Sw_s) * (1 - Sro - Sw_s) * mult_o 
		   + (Sw_s - Srw) * (Sw_s - Srw) * mult_w ) * dx_dy;
      aP (i, j) = aE (i, j) + aW (i, j) + aN (i, j) + aS (i, j);
    }
    applyBoundaryConditions2D();
    return 0;
}


} // Tuna namespace


#endif //_BLIP2_H_

















